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1.0  Introduction 


The  problem  of  constructing  interpolation  functions  for  sets  of 
scattered  data  in  two  independent  variables  has  been  treated  in  many  papers. 
The  survey  paper  on  approximations  to  multivariate  functions  by  Schuoaker  [20] 
contains  extensive  references.  Other  survey  papers  are  by  Barnhill  [3] 
and  Sabin  [19].  The  present  author  has  surveyed  and  tested  a  number  of 
algorithms  for  solution  of  the  problem  [10],  [11].  Conceptually  the 


problem  is  quite  simple:  Given  points  (x^.y^.f^),  k  *  1,  . ..,  N,  with 
distinct  ) »  construct  a  function  F(x,y)  so  that  FCx^.y^)  =*  f^, 

k  *  1,  ...»  N.  Generally  one  wants  to  have  a  smooth  interpolant,  F(x,y), 


in  the  sense  that  low  order  partial  derivatives  are  everywhere  continuous. 


This  is  complicated  for  large  sets  of  data  by  the  fact  that  the  interpolant 

(in  a  practical  sense,  to  be  computable)  must  be  local,  so  that  its  value 

at  some  point  (x,y)  depends  only  on  values  for  which  (x^,y^)  is 

"close"  to  (x,y) .  A  general  framework  for  a  class  of  such  methods  is 

given  in  [9],  and  we  will  discuss  it  briefly  in  Section  2. 

While  a  large  number  of  ideas  have  been  proposed  for  solution  of 

the  problem,  a  much  smaller  number  of  working  computer  programs  are  readily 

available.  These  include:  (1)  a  method  based  on  finite  element  functions 

defined  over  triangles,  due  to  Akima  [1],  [2],  a  version  of  which  is 

* 

available  in  the  IMSL  library  under  the  name  IQHSCV,  (2)  a  program  based 
on  a  similar  idea,  due  to  Lawson  [14],  (3)  two  programs  based  on  weighted 
local  approximations  by  quadratic  functions,  due  to  Franke  and  Nielson  [12]. 
A  program  by  Little  [15],  and  another  by  Nielson  [18],  both  based  on  finite 


* 

International  Mathematical  and  Statistical  Libraries,  7500  Bellaire  Blvd., 
Houston,  TX  77036. 
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element  functions  defined  over  triangles,  will  probably  be  available 
by  press  time.  All  of  these  programs  have  been  tested  by  the  author 
[10],  [11],  and  most  perform  adequately  in  a  variety  of  cases;  None  of 
them  seems  to  have  a  clear  edge  over  all  the  others,  or  to  be  entirely 
satisfactory.  For  certain  applications,  each  has  its  good  points.  The 
choice  of  a  method  for  most  users  will  be  based  on  subjective  criteria 
which  vary  from  person  to  person  and  application  to  application.  It  is 
not  surprising  this  is  the  case  in  two  variables  since  it  is  also  the  case 
in  one  variable,  although  perhaps  to  a  lesser  extent. 

The  purpose  of  this  paper  is  to  document  an  alternative  scheme  which 
performed  comparably  well  in  the  previously  mentioned  tests.  In  this  way 
the  computer  program  will  be  brought  to  the  attention  of  potential  users, 
who  may  request  it  from  the  author.  The  method  will  supplement  the  current¬ 
ly  available  codes  in  that  it  is  based  on  a  different  approach.  It  is 
anticipated  that  it  will  find  application  and  approval  in  a  variety  of  areas. 

The  theoretical  background  for  the  method  is  discussed  in  Section  2, 
while  details  of  the  program  are  outlined  in  Section  3.  Section  4  gives 
information  concerning  usage  of  the  program.  Several  examples  are  given 
in  Section  5,  including  perspective  plots  of  some  surfaces  generated  by  the 
program.  Included  is  an  example  of  how  the  variation  of  a  parameter  in 
the  method  affects  the  surface.  General  guidelines  for  choice  of  this 
parameter  are  given,  even  though  the  suggested  value  usually  leads  to 
satisfactory  results.  A  general  discussion  of  the  features  of  the  method 
is  given,  along  with  general  guidance  for  use  of  the  program. 

2.0  Theoretical  Background 

The  general  idea  encompassing  this  scheme  and  many  others  are  given 
in  [9],  Consider  that  the  points  (x^.y^.f^) ,  k  *  1,  ...,  N  are  given. 
Briefly,  local  approximations  to  the  data  are  constructed,  and  these  are 
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then  blended  together  by  using  weighting  functions  which  form  a  partition 
2 

of  unity  on  R  .  We  pause  to  give  the  necessary  notation  and  definitions. 

A  set  of  functions,  w^(x,y),  i  *  1,  . . . ,  m  is  said  to  form  a  partition 

m 

of  unity  if  each  w.(x,y)  £  0  and  £  w. (x,y)  =  1.  The  w  will  be  called 

i-1  1 

weight  functions.  Let  the  support  of  w^  be  denoted  by  -  closure  {(x,y): 

w^x.y)  /  0}.  Let  »  {k:  (x^,y^)  £  Interior  (S^)}.  Now  suppose  that  the 

functions  Q^Cx.y),  i  *  1,  . . . ,  m,  are  defined  on  and  have  the  property 

that  they  interpolate  the  data  whose  (x,y)  coordinates  are  in  Interior  (S4), 

i.e.,  if  k  £  1^,  then  (^(m^.y^)  *  f^.  These  functions  0^  will  be  called 

m 

local  approximations.  We  then  consider  the  function  F(x,y)  *  J  w  (x,y)Q. (x,y) . 

i-1  1 

Its  properties  are  summarized  in  the  following. 

m 

Proposition.  The  function  F(x,y)  -  £  w  (x,y)0. (x,y)  has  the  following 

i-1 

properties: 

(1)  Interpolation;  F(xk>yk)  -  ffe,  k  -  1,  ....  N. 

(2)  Smoothness;  F(x,y)  is  at  least  as  smooth  as  the  wi  and  Q±,  e.g. , 
if  all  of  the  functions  v^,  Qi>  i-1,  ...,  m  have  continuous 
first  derivatives,  so  does  F(x,y). 

(3)  Local  dependence  on  the  data;  Let  Cx,y)  be  fixed,  and  let 

Jx  y.  -  (i:  w^(x,y)  +  0},  then  F(x,y)  depends  only  on  the  (x^jy^.f^) 

points  for  which  k  £(  u  I. )  j { i :  some  0.,  j  e  J  depends  on 

ieJ  1  3 

x.y 

^Xi’yi,fi^''  In  ParCicular>  we  have  f(x,y)  -  \  Wi(x,y)Qi(x,y) . 

^  icJ 

x.y 

These  properties  are  essentially  observations,  but  form  the  basis  for 
construction  of  appropriate  weight  functions  which  will  allow  easy  determin¬ 
ation  of  the  set  J  .  Our  construction  vields  a  set  of  at  most  four  nonzero 
x,y 

terms  in  the  sum  defining  F(x,y).  This  provides  a  considerably  faster  process 
during  the  evaluation  of  the  interpolant  than  was  possible  in  the  choices 


previously  considered  [9]. 


A 


r 


It  has  beer  Implied,  but  is  not  necessary  for  the  proposition  to 

hold,  that  many  of  the  weight  functions,  w^  have  finite  support.  In  order 

for  the  local  approximations  Q^(x,y)  to  actually  be  local,  this  will 

likely  be  the  case.  Therefore  we  think  of  weight  functions  whose  support, 

S^,  contains  relatively  few  Cx^.y^)  points,  and  whose  support  is  often 

2 

local.  In  order  for  the  interpolant  to  be  defined  on  all  of  R  ,  some  sets 
must  not  be  finite,  and  typically  would  contain  points  (x^.y^)  on  or 
near  the  boundary  of  the  convex  hull  of  the  set  of  points  {(x^,y^)}.  With 
this  framework  and  ideas  in  mind,  we  are  ready  to  discuss  the  specific 
details  of  the  algorithm. 


3.0  Details  of  the  Algorithm 


3.1  Choice  of  Weight  Functions 

While  the  choice  of  weight  functions  was  allowed  to  determine  the 
support  regions  in  the  discussion  of  the  previous  section,  it  is  more 
convenient  to  proceed  from  support  regions  to  weight  functions  in  the 
actual  application.  The  use  of  support  regions  which  are  rectangles  have 
specific  advantages  in  terms  of  controlling  the  number  of  support  regions  in 
which  a  particular  point  (x,y)  lies,  as  well  as  simplifying  determination 
of  those  regions. 

Let  n  and  n  be  given  positive  integers,  and  let  finite  values  of 
x  y 

*0  <  *1  <  *2  <  "  Xn  <  xn  +1  and  y0  <  yl  '  y2  "  **•  "  ?n  <  yn  +1  be 

xx  y  y 


given.  For  each  i  ■  1,  2,  ....  n  and  j  ■  1,  2,  ...,  n  let  r,  .  * 

x  y 

[xi-i*  Xi+1]  *  [yj-i’  Vi]- 

2  3 

Let  Hj(s)  *  1  -  3s  +  2s  ,  the  Hermite  cubic  satisfying  H^CO)  ■  1, 
Hj(l)  •  H^(0)  ■  *  0.  We  then  define  piecewise  cubics  with  continuous 


first  derivatives,  which  are  nonzero  only  on  two  adjacent  intervals,  and 
satisfy 


The 


v.(x.)  -  5t.  ,  i.J  -  1.  2,  ....  ^ 


uj(yi>  "  5ji 


i,j  -  1,  2,  ....  n 


In  particular  we  take 
'  1 


vl(x)  ■  < 


h3| 

^  0 


x  -  x. 


x2  -  xx 


X  <  X, 


X^  s  X  <  X, 


X  2  X„ 


vi(x) 


1  -  Vj.^x) 


X  -  X, 


0 


xi+l  '  Xi 


V. 


X  <  X 


i-1 


xi-l  “  X  xi 


xt  .  X  <  x1+1 


X  i  X 


i+1 


for  1=2,  ....  n^  -  1 


r 


X  <  X 


n  -1 
x 


vn  (x)  »  -S 
x 


1  -  -1»  •  *n  -1  S  *  ‘  Xn 

X  X 

1  .  X  >  x_ 


Accer-Flon  "or 
NTIS  P'm&i 
DTIC  5MB 
Unannounced 
Juttif  lent  ion. . 

By -  - 

Di stribrt 1  on/ 

Availnl'.  l  j  t.v 

a-.-.; i  V.- . 

[D 1 .. t  :*p.  j  •  ; 


u^(y)  are  dual. 


Then  if  we  define 


WL  ^(x.y)  »  vi(x)ui(y) 


i-1,  ...,nx,  j  -  1.  •••»  n 
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it  is  easilv  observed  that  w.  .  has  support  r.  ,,  except  for  when  i  *  1  or  n  , 

3-.J  i,j  x’ 

or  j  *  1  or  n^,  when  the  support  extends  to  30  in  one  or  both  variables.  We 

denote  the  support  of  W.  ,  bv  R.  ..  Further,  we  note  that  the  W,  .  form  a 

i.j  ‘  i.J  i.J 

o 

partition  of  unity  on  R“. 

Since  the  x^  and  v^  values  give  rise  to  a  grid  on  the  plane,  we  call 

them  grid  values.  The  choice  of  grid  values  x^  and  y  depend  on  the  data. 

They  may  be  specified  by  the  user,  but  the  default  option  is  for  them  to  be 

determined  by  the  program.  In  the  default  mode  the  user  specifies  a  parameter, 

NPPR  (for  number  of  points  per  region),  the  suggested  value  being  about  10. 

The  program  will  then  determine  the  grid  values  so  that  the  anticipated 

number  of  (x,  ,y,  )  points  in  each  rectangle  R  will  be  approximated  NPPR. 

x  k  i ,  j 

For  data  which  Is  not  somewhat  uniformly  distributed  the  actual  numbers  may  vary 
considerably.  However,  for  most  situations  we  have  encountered,  the  process 
is  quite  adequate. 

An  equal  number  of  grid  values  is  taken  in  each  direction,  i.e.,  n^  *  n^. 

Because  we  want  NPPR  points  per  rectangle,  each  subrectangle  (x^.x^^)  x  (y^.y^^) 

1  ">  1 
should  have  NPPR  points.  Since  n^  ■  n^,  we  want  (n^  +  1)“  —  NPPR  »  N, 

1/2 

leading  us  to  take  n^  to  be  the  nearest  integer  to  (4N/NPPR)  -  1. 

The  grid  values  x^,  are  now  chosen  so  that  there  are  approximately  equal 
numbers  of  points  from  the  values  x^,  k  ■  1,  . . . ,  N  in  each  interval 
(Xf,  x^+1) .  Specifically,  let  x^  denote  the  x^  values  arranged  in  non¬ 
decreasing  order.  Consider  the  points  (0,x^),  (l.Xj),  ....  (N-l,  x^) ,  and 
let  g(t)  be  the  piecewise  linear  interpolant  for  them.  Divide  the  interval 


N-l 

(0,N-1)  into  n  +  1  3ubintervals  of  length  A  ■  — -r- 
x  n  +1 

x 


The  values  of  x^  are 


determined  by  taking  them  to  be  the  values  of  g  at  the  endpoints  of  the 
subintervals,  i.e.,  x^  *  g(i  ),  i  *  0,  1,  2,  ...,  n^+  1.  The  y  are  determined 
in  dual  fashion.  This  process  results  in  the  grid  values  and  hence  the 


rectangles,  being  symmetric  if  the  (x^.y^)  points  are  symmetric.  In  addition, 
the  relative  positions  of  grid  values  are  unchanged  by  linear  displacements 
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and  stretching  in  each  variable. 


When  chosen  in  the  above  fashion,  the  location  of  the  lines  is  not 
a  local  process  in  the  sense  that  insertion  of  an  additional  point  will 
change  the  boundaries  of  all  of  the  rectangles.  While  one  could  argue  that 
then  the  scheme  is  not  local,  we  take  the  view  that  the  idea  of  local 
determination  of  the  interpolant  is  most  important  in  the  evaluation  phase. 

The  determination  of  parameters  in  the  scheme  (here,  the  x^  and  v^)  may  be 
a  global  process.  Of  course,  if  the  user  specifies  the  grid  values,  he  will 
likely  be  using  a  global  process  co  choose  them. 

3.2  Choice  of  Local  Approximations 

The  only  constraint  on  the  local  approximations  is  that  they  interpolate 
the  appropriate  points,  and  that  they  have  continuous  first  derivatives 
(at  least)  to  assure  a  smooth  interpolant.  In  the  previously  mentioned  tests 
conducted  by  the  author,  a  number  of  global  interpolation  schemes  for  scattered 
data  were  considered.  In  principle,  any  of  these  might  be  used.  The  choice 
here  was  made  for  two  reasons,  (1)  the  method  scored  very  well  in  the  tests, 
and  (2)  the  method  has  an  elegant  and  well  developed  mathematical  theory  which 
also  has  direct  application  to  some  engineering  problems. 

The  local  approximations  used  in  this  algorithm  are  the  thin  plate  splines 
first  mentioned  by  Harder  and  Desmairis  [13],  with  theoretical  developments 
by  Duchon  [4-7]  and  Meinguet  [16],  [17].  It  was  first  developed  as  the 
solution  to  the  problem  of  a  thin  plate  which  is  forced  to  pass  through 
certain  points  (the  interpolation  points)  by  application  of  point  loads.  For 
our  purposes  it  is  sufficient  to  know  that  the  approximation  is  of  the  form 


where  I  • 


(y  -  yk>2 


r  2 

Q(x,y)  =*  I  A^d^  log  dk  +  a  *  bx  +  cv  , 
kel 

2  2 

(k:  Q  is  to  take  on  the  value  f^  at  and  dk  *  (x  -  x^)-  'u 

The  coefficients  A.  ,  and  a,  b,  and  c  are  determined  by  the  linear 


♦ 


{ 
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system  of  aquations 


n  z 

)  A^d^  log  d,  +  a  +  bx  +  cv 


Cx.y) 


(xl’Vi) 


fi  *  ieI 


l  \  m  0 

kel 

l  Vk  * 0 

kel  w  K 

.  L  Vi  ■ 0 


The  geometric  effect  of  the  last  three  equations  is  to  suppress  all  terns 
in  the  approximation  which  grow  faster  than  linear  as  distance  from  the 
interpolation  points  is  increased.  A  linear  system  of  order  equal  to  the 
number  of  interpolation  points  plus  three  must  be  solved.  To  be  nonsingular 
there  must  be  at  least  three  noncollinear  points  among  the  (x^,y^),  kel. 

The  system  is  symmetric,  but  not  positive  definite.  While  an  equation 
solver  designed  for  such  systems  could  be  used,  we  have  found  a  general 
purpose  solver,  the  DECOMP/SOLVE  subroutines  of  Forsythe,  Malcolm,  and  Moler  [3] 
has  given  more  reliable  results. 

It  is  easily  observed  that  the  local  interpolant  has  continuous  derivatives 
of  all  orders  except  at  the  data  points,  (x^.y,^),  where  a  logarithmic 
singularity  occurs  in  the  second  derivatives.  The  interpolant  has  linear 
precision,  that  is,  if  the  (x^.y^.f^)  points  all  lie  on  a  linear  function,  the 
interpolant  will  reproduce  it. 

While  the  thin  plate  spline  is  invariant  with  respect  to  scaling,  trans¬ 
lation,  and  rotation  (not  all  of  this  is  obvious),  the  condition  number  of 
the  coefficient  matrix  for  the  system  of  equations  is  dependent  on  scaling. 

To  minimize  difficulties  with  that,  and  to  remove  the  effect  of  scaling  the 

variables  by  different  amounts,  we  transform  each  rectangle  r  .  onto  the  unit 

♦  J 

2 

square  [0,1]“,  before  the  local  approximation  0  is  computed. 

*  J 


s 


Ic  remains  to  specifv  the  process  for  deterwi nine  the  points  (x,  .  , f,  ) 

ic  tc  k 

co  be  Interpolated  by  the  thin  place  spline  local  approximation.  Experience 
has  shown  that  it  is  advantageous  to  include  more  points  than  is  necessarv, 
i.e.,  (x,  ,v,  )  which  are  outside  cf  R,  ..  This  tends  to  yield  a  better  dans- 

rC  1  ,  J 

ition  between  local  approximations  than  when  only  necessary  points  are  included. 

a 

The  set  of  (x^y^)  points  transformed  into  the  rectangle  R  =  [-.1125,  1.1125]“ 
by  the  transformation  taking  d  ^  onto  [0,1]“  are  included.  Let 


I.  .  =  (k: 

i.l 


*k  '  Xi-1 


y,  -  v 

k  ->-1 


x..t  -  x.  ,  v.,,  -  v. 

i+l  i-l  -j+1  -  j-1 


£  R; . 


This  gives  the  basic  set  of  interpolation  points  for  the  local  approximation 

Q.  .  associated  with  the  rectangle  R.  . .  Undar  certain  conditions  there  mav 
i,3  i,j 

be  fewer  than  the  necessarv  three  indices  in  I.  ..  When  this  hapoens,  the 

i.J 

set  I.  .is  augmented  by  including  as  interpolation  points  the  necesary 
1 1 J 

number  of  closest  points  to  the  rectangle  S.  ^  (in  the  norm) ,  after  the 
points  (x^.y^)  have  undergone  the  transformation  to  the  unit  square.  The 
minimum  number  of  points  per  rectangle  is  a  variable,  MINPTS .  This  has  been 
sec  to  3,  but  may  be  increased  if  it  seems  desirable. 

After  the  interpolation  points  for  each  local  approximation  have  been 
determined,  the  local  thin  plate  splines,  can  be  determined  bv 

calculating  the  coefficients.  This  yields 

Q,  .(x,v)  =  /  A.  .  ,  d,'2  log  d,’  +  a.  .  +  b.  .  x'  +  c.  .  v' 

l.J.K  k  ^  l.J  1,1  • 

where  the  primes  denote  coordinates  and  distance  after  the  transformation  of 

d  j  to  the  unit  square 

3.3  Properties  of  the  Interoolant 


The  overall  interpolant  is  of  the  form 
nx  n 

F(x,y)  -  l  f  W  (x,v)Q  (x.y), 
i-l  j-1  1>;1 
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and  as  noced  previously,  there  are  at  most  four  nonzero  terms  in  the  sum. 

Which  terms  are  nonzero  is  easilv  determined.  If  x  2  x  ,  set  i'  ■  n  ; 

n  x 

x 

Otherwise  let  i'  be  the  smallest  index  so  that  x  ,  >  x.  Determine  j  '  in 

dual  fashion  from  v  and  the  v,'s.  Then,  the  four  nonzero  terms  have  indices 

'  J 

(i’.j'Mi'+l.J'),  (i'.j’+l),  and  (i'+l.j'+l).  If  i'  =  0  or  i'  *  nx>  the 
terms  involving  i'  or  i'+l,  respectively,  do  not  appear,  and  similarly  if 
j'  =  0  or  j ’  ■  n  ,  the  terms  involving  j*  or  j'+l,  resoectively ,  do  not 
appear. 

In  dition  to  the  properties  outlined  in  Section  2,  certain  other 
properties  hold.  The  approximation  is  invariant  under  translation  and 
stretching  (independently  in  each  variable).  It  has  symmetry  with  respect 
to  planes  parallel  to  coordinate  planes  whenever  the  data  has  that  symmetry. 
The  approximation  is  not  invariant  under  rotations,  however,  since  the 
rectangles  depend  on  the  individual  coordinates  of  the  data  points. 

The  approximation  has  continuous  first  derivatives,  and  -jump  dis¬ 
continuities  in  the  second  derivatives  across  grid  lines,  as  well  as 
logarithmic  singularities  in  the  second  derivatives  at  the  data  points. 

Plots  of  the  surfaces  generally  appear  to  be  quite  smooth,  however.  Since 
the  local  approximations  have  linear  precision,  the  overall  approximation 
also  has  linear  precision,  i.e.,  if  the  data  lies  on  a  linear  function,  the 
interpolant  is  a  linear  function. 

4.0  The  Computer  Program 

The  overall  hierarchy  of  the  subroutines  is  given  in  Diagram  1,  which 
shows  the  communication  links  between  them.  No  COMMON  is  used  as  array 
sizes  are  problem  dependent  and  thus  specified  by  the  user.  We  brieflv 
discuss  them,  mentioning  important  parameters. 
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Diagram  1. 


Calling  Program 

This  program  is\ supplied  by  Che  user  and  must  supply  the  (x.^  ,y  ,  f^) 

points,  plus  two  arrays  of  points  xO . ,  and  yO .  for  the  grid  of  points 

i  1  J 

(xO^.yOj)  at  which  the  interpolant  is  to  be  evaluated.  In  addition,  the 

user  must  supply  two  workspace  arrays,  IWK  and  WK  in  which  information 

calculated  during  preprocessing  (e.g.,  I  .,  x.,  y  ,  and  coefficients  for 

i  >  3  1  J 

the  local  approximations),  is  stored,  and  an  arrav  FO  for  the  returned 
interpolant  values. 

The  amount  of  storage  required  for  the  arrays  IWK  and  WK  is  not 
known  a  priori.  The  estimated  space  required  is  about  6N  for  IWK  and 
about  7N  for  WK.  Table  1  gives  exact  results  for  several  different  sized 
problems  based  on  random  (x,y)  points.  Oddly  distributed  point  sets  may 
result  in  somewhat  more  storage  being  required.  In  any  case,  the  precise 
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number  of  locations  required  is  returned  to  the  calling  program  from 


Subroutine  LOTPS,  the  only  routine  referenced  by  the  user.  If  an  in¬ 
sufficient  number  are  allowed,  an  error  return  occurs. 

Under  the  usual  option,  the  user  specifies  NPPR  >  0,  the  suggested 
value  being  10.  If  the  user  wishes  to  specify  the  grid  lines,  he  may 
do  so  by  setting  JJPPR  ■  0,  and  then  giving  grid  line  information  in  the 
arrays  IWK  and  WK,  as  explained  in  the  argument  description  for  Subroutine 

LOTPS.  Typically  one  should  take  Xq  *  m^Ln  x,^,  x^  »  a^x  x^ ,  and  the  dual 

x 

in  y.  This  is  not  necessary,  although  all  points  to  be  interpolated  should 

lie  in  [x.,  x  . .  J  [y.,  v  To  prevent  different  scaling  (intemallv) 

u  n  ■ri  U  '  n 

x  y 

in  the  two  variables,  a  square  grid  covering  the  (x^,/,^)  points  could  be 
specified. 

Subroutine  LOTPS 

This  subroutine  provides  the  interface  between  the  user's  program  and 
the  sec  of  routines  implementing  the  method.  Generally,  LOTPS  sets  up 
storage  areas  In  the  arrays  IWK  and  WK,  determines  parameters  required  by 
other  subroutines,  and  calls  other  subroutines  to  (1)  generate  the  grid.  If 
necessary,  (2)  determine  the  interpolation  points  for  the  local  approximations, 
(3)  compute  coefficients  for  the  local  approximations,  (4)  evaluate  the 
interpolant  at  a  grid  of  points. 

Subroutine  GRID 

This  subroutine  selects  values  of  x^  and  v^  in  accordance  with  the 
discussion  in  Section  3. 

Subroutine  LOCLIP 


This  subroutine  determines  the  interpolation  points  for  each  local 


interpolant,  in  accordance  with  Section  3. 

Subroutine  CLOTPS 

This  subroutine  generates  the  system  of  equations  for  the  coefficients 
of  the  local  approximations  and  calls  an  equation  solver  to  obtain  them. 


Internally  the  routine  has  a  maximum  of  30  local  interpolation  points. 

This  can  easily  be  altered  by  changing  two  statements,  as  noted  in  the 
program  listing. 

Subroutine  EVLTPS 

This  subroutine  evaluates  the  interpolant  on  the  grid  of  points  specified 
by  the  user.  Use  of  a  grid,  when  that  is  required,  facilitates  the  process 
of  locating  the  rectangle  in  which  a  particular  evaluation  point  is  located. 
The  xO^  and  yO^  values  should  be  in  increasing  order  for  maximum  efficiency. 
Evaluation  time  for  a  grid  of  points  should  be  nearly  independent  of  N  for 
large  N,  which  is  borne  out  by  the  timing  information  in  Table  2. 


5.0  Examples  and  Observations 

Example  1.  This  example  shows  a  typical  local  approximation  function  and 
is  for  the  data  in  Table  1.  This  function  is  a  "cardinal"  function  for  the 
first  point,  and  as  such  shows  the  effect  of  a  nonzero  value  at  a  single  point 


on  the  interpolant.  The  plotted  surface,  shown  in  Figure  1,  is  over  [0,1]“. 

Example  2.  This  example  is  given  to  show  a  surface  generated  from  100 

randomly  generated  points  with  the  function  value  being  obtained  from  an 

explicitly  given  function.  The  surface  was  generated  using  NPPR  *  10,  has 

rms  error  of  about  .3%,  and  is  virtually  indistinguishable  from  a  plot  of  the 

parent  surface.  It  is  shown  in  Figure  2. 

Examples  3 ,  4 ,  5 .  These  examples  use  the  same  set  of  60  points  lying 
1  19  2 

in  the  square  (  -  ’  ~i8  ^  ’  an<*  c^osen  a  pseudorandom  number  generator. 

The  function  is  explicitly  given  by  f(x,y)  *  .1  +  ,  and  is 

(x+y ) 

shown  in  Figure  3.  Figures  4,  5,  and  6  show  the  interpolant  over  the  square 
[0,1]^  for  NPPR  values  of  6,  10,  and  15. 

From  these  examples  we  see  that  NPPR  *  10  works  well,  not  too  much 
difference  i3  observed  when  NPPR  is  increased,  but  NPPR  *  6  gives  a  less 
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smooch  appearing  surface.  The  smaller  MPPR  is,  the  more  localized  the 

surface  becomes  (although  MPPR  ■  1  will  probably  be  the  least  value 

practicable,  and  some  local  interpolants  will  likely  become  planes  due  to 

the  minimum  of  three  interpolation  points  being  reached).  In  line  with 

this  comment,  very  smooth  surfaces  with  small  gradients  will  probably  be 

^amenable  to  larger  MPPR,  while  surfaces  with  large  gradients  may  be  best 

approximated  by  taking  MPPR  smaller,  thus  localizing  the  behavior. 

Table  2  gives  the  results  of  a  series  of  problems  with  various  numbers 

of  points  and  values  of  MPPR.  The  data  points  were  chosen  bv  a  pseudorandom 

1  19  2 

number  generator  in  the  square  [  -  yg  •  jg-  ]  ,  and  the  approximation  was 


evaluated  on  a  33x33  grid  (of  1089  points  total)  on  [0,1  P.  We  observe  that 
increasing  MPPR  Increases  execution  time  while  decreasing  the  amount  of 
storage  needed.  Preprocessing  time  should  be  about  proportional  to  >T ,  but 
apparently  there  is  a  strong  linear  component  for  small  M.  Preprocessing 
time  should  also  be  about  proportional  to  MPPR,  although  this  is  not  readilv 
apparent.  Evaluation  time  should  be  nearly  independent  of  M  for  large  S,  and 
proportional  to  MPPR,  which  is  approximately  shown. 

The  automatic  grid  selection  process  works  well  when  the  data  is 
fairly  uniform  in  (x,v)  and  lies  nearly  in  a  square  region.  If  the  data  is 
very  irregular,  or  lies  in  an  oblong  rectangular  area,  it  will  probably  be 
useful  to  explore  results  with  a  user  specified  grid  to  obtain  better  coverage 
by  rectangles  which  are  not  too  oblong.  Limited  experience  has  been 
accumulated  in  these  situations. 


6.0  How  to  Obtain  this  Program 

A  copy  of  this  program,  written  in  Fortran,  and  including  a  sample 
driver  program,  can  be  obtained  from  the  author.  To  do  so,  send  a  (short) 
1/2  inch  tape  to  the  author  indicating  the  format  desired. 
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X 


y 


f 


0.35 

0.35 

0.5 

-0.05 

0.25 

0.0 

0.10 

-0.05 

0.0 

0.50 

0.05 

0.0 

0.00 

0.90 

0.0 

0.30 

0.70 

0.0 

0.60 

0.50 

0.0 

0.90 

0.00 

0.0 

0.40 

1.05 

0.0 

0.85 

0.80 

0.0 

1.05 

0.20 

0.0 

1.10 

1.10 

0.0 

Table  1: 

Data  for  "Cardinal" 

function. 
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>n»PR 

N 

n  *  n 

X  V 

NWKU 

NIWKU 

Time 

(Preprocessing) 

Time 

(Evaluation) 

6 

60 

5 

334 

273 

4.2 

14.3 

6 

100 

7 

619 

506 

7.5 

17.6 

6 

500 

18 

3596 

2893 

70.5 

21.2 

6 

1000 

25 

7441 

6140 

201.5 

21.3 

10 

60 

4 

284 

243 

5.3 

18.0 

10 

100 

5 

488 

427 

9.9 

26.4 

10 

500 

13 

3014 

2649 

73.3 

29.8 

10 

1000 

19 

6565 

5804 

202.6 

31.7 

15 

60 

3 

224 

199 

6.8 

23.0 

15 

100 

4 

414 

373 

12.4 

32.8 

15 

500 

11 

2856 

2591 

91.6 

39.3 

15 

1000 

15 

5920 

5439 

258.5 

40.6 

Table  2:  Storage/Times  for  Various  Size  Problems. 
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Figure  1:  Cardinal  ~>"ction 


Figure  2:  Saddle  Function 
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SUBROUTINE  L.OTPS  (MODE. NPPR# NPI » XI , YI » FI , NXO. XO, NYO. YO, I WK . HI WK,  LOT 
NIMKU/ WK. NWK. NWKU. FO. KER  )  LOT 
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